
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

#include "util.h"
#include "ibm2d.h"
#include "nsmg.h"
#include "boundary.h"
#include "dlist.h"
#include "output.h"

#define REAL double


void print_data_uvp(REAL **u, REAL **v, REAL **p, int count, char *buf)
{

	char bufferu[20], bufferv[20], bufferp[20];
	int i, j;

	FILE *fu, *fv, *fp;

	if (count < 10) {
		sprintf(bufferu, "./%s/u000%d.m", buf, count);
		sprintf(bufferv, "./%s/v000%d.m", buf, count);
		sprintf(bufferp, "./%s/p000%d.m", buf, count);
	}
	else {
		sprintf(bufferu, "./%s/u00%d.m", buf, count);
		sprintf(bufferv, "./%s/v00%d.m", buf, count);
		sprintf(bufferp, "./%s/p00%d.m", buf, count);
	}

	fu = fopen(bufferu, "w");

	iloop{
		jloop {
			fprintf(fu, "  %16.14f", 0.5*(u[i - 1][j] + u[i][j]));
		}
		fprintf(fu, "\n");
	}
	fclose(fu);

	fv = fopen(bufferv, "w");
	iloop{
		jloop {
			fprintf(fv, "  %16.14f", 0.5*(v[i][j - 1] + v[i][j]));
		}
		fprintf(fv, "\n");
	}
	fclose(fv);

	fp = fopen(bufferp, "w");
	iloop{
		jloop {
			fprintf(fp, "  %16.14f", p[i][j]);
		}
		fprintf(fp, "\n");
	}
	fclose(fp);

	return;
}


void print_data_bdpo_normal_curv(REAL **p, REAL **normal, REAL *curvature, int nt, int count, char *buf_folder)
{
	int k;
	char buf[20];
	FILE *fp;

	if (count < 10)
		sprintf(buf, "./%s/bdpo000%d.m", buf_folder, count);
	else
		sprintf(buf, "./%s/bdpo00%d.m", buf_folder, count);

	fp = fopen(buf, "w");

	for (k = 0; k <= nt; k++) {
		fprintf(fp, "%f %f %f %f %f \n", p[k][1], p[k][2], normal[k][1], normal[k][2], curvature[k]);
	}

	fclose(fp);

	return;
}


void print_data_LB(REAL **p, int Np, REAL time, char *buf)
{
	extern REAL pi;

	int i, j, index_max_d;
	REAL center[3], *d, max_d, min_d, theta, D, L, B;
	FILE *my;

	d = dvector(1, Np);

	make_center(p, Np, center);

	for (i = 1; i <= Np; i++) {
		d[i] = 0.0;
		for (j = 1; j <= 2; j++) {
			d[i] += pow(p[i][j] - center[j], 2);
		}
		d[i] = sqrt(d[i]);
	}

	max_d = min_d = d[1];
	index_max_d = 1;
	for (i = 2; i <= Np; i++) {
		if (d[i] > max_d) {
			max_d = d[i];
			index_max_d = i;
		}
		if (d[i] < min_d)
			min_d = d[i];
	}

	L = max_d;
	B = min_d;
	D = (L - B) / (L + B);
	theta = atan2(p[index_max_d][2] - center[2], p[index_max_d][1] - center[1])*180.0 / pi;

	my = fopen(buf, "a");
	fprintf(my, "%e %e %e %e %e \n", time, theta, D, L, B);
	fclose(my);

	free_dvector(d, 1, Np);
}
